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Abstract 

Uncertainty quantification (UQ) techniques are frequently used to 
ascertain output variability in systems with parametric uncertainty. 
Traditional algorithms for UQ are either system-agnostic and slow 
(such as Monte Carlo) or fast with stringent assumptions on smooth- 
ness (such as polynomial chaos and Quasi- Monte Carlo). In this work, 
we develop a fast UQ approach for hybrid dynamical systems by ex- 
tending the polynomial chaos methodology to these systems. To cap- 
ture discontinuities, we use a wavelet-based Wiener-Haar expansion. 
We develop a boundary layer approach to propagate uncertainty through 
separable reset conditions. We also introduce a transport theory based 
approach for propagating uncertainty through hybrid dynamical sys- 
tems. Here the expansion yields a set of hyperbolic equations that are 
solved by integrating along characteristics. The solution of the par- 
tial differential equation along the characteristics allows one to quan- 
tify uncertainty in hybrid or switching dynamical systems. The above 
methods are demonstrated on example problems. 

1 Introduction 

Uncertainty Quantification (UQ) is an area of mathematics that is used to 
quantify output distributions given parametric uncertainty. Traditional ap- 
proaches include Monte Carlo and Quasi-Monte Carlo methods [1], response 
surface methods [21 13] as well as polynomial chaos and probabilistic colloca- 
tion based approaches [Ij. The polynomial chaos approach for uncertainty 
quantification was originally proposed by Norbert Wiener [5]. Assuming 
that one is given input uncertainty in the form of distributions associated 



with various parameters of the system, polynomial chaos/probabilistic col- 
location methods provide an approach for fast uncertainty quantification 
under the assumption of smooth dynamics. In particular, polynomial chaos 
provides exponential convergence for smooth systems and processes with 
finite variance [1]. Polynomial chaos based methods have been used for a 
multitude of apphcations, see [6l [71 El El [13 III Ell 13] for examples. Note 
that, depending on the application, one can combine various UQ approaches. 
For example, a combination of polynomial chaos and the response surface 
methodology has been used to develop probabilistic collocation methods for 
discrete distributions in [TO] . 

In this work, we focus on developing UQ techniques for hybrid dynamical 
systems. Hybrid dynamical systems theory is used to model systems with 
both discrete and continuous dynamics [H] . Examples include the bouncing 
ball automaton [15], biological networks [16\ 117] , air traffic management 
systems [18], communication networks [19], elevators, and robotics, to name 
a few. These systems frequently display rich dynamics not seen in continuous 
systems. For example, Zeno behavior in hybrid systems is characterized by 
an infinite number of discrete switches in finite time [151 [20] . Hybrid systems 
can be particularly challenging from an analysis standpoint since traditional 
techniques, such as polynomial chaos based methods, assume smoothness, 
rendering them inapplicable. 

In this work, we develop polynomial chaos and transport theory based 
methods for propagating uncertainty through hybrid systems. We assume 
that the domains associated with different modes of operation of the hybrid 
system do not overlap. We demonstrate that, by integrating over appropri- 
ate time-varying regions, one can extend the polynomial chaos framework 
to hybrid dynamical systems. We resolve the issue of state resets jl4] in 
the separable case by using boundary layer approximations. To capture 
the discontinuities in the probability distributions of the output variables, 
we use a Haar- wavelet expansion [2T]. This expansion has previously been 
used in the polynomial chaos setting to propagate uncertainty through dy- 
namical systems close to bifurcation points [H \22[ 123] . Here we develop a 
methodology to propagate uncertainty through systems with discontinuities 
in dynamics and output along with state resets. We also develop a trans- 
port theory based approach that allows one to propagate the uncertainties 
through the various modes of the hybrid dynamical system. 

Our paper is organized as follows: in section [2] we define hybrid dynam- 
ical systems and the problem of uncertainty quantification. In section [3] we 
first construct the framework for polynomial chaos in the hybrid dynamical 
system setting (13. ID . We then demonstrate the Haar wavelet expansion for 
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hybrid polynomial chaos in 13.21 The handling of state resets is considered 
in section [3^ Finally, the results on hybrid polynomial chaos are presented 
in 13.41 The transport operator theory based method for propagating un- 
certainty through hybrid dynamical systems is developed in section [J] and 
conclusions are drawn in section O 

2 Problem definition 

Let S = {q, X, /, x(0), D, E, G, R) denote a hybrid system 5, where 



R: E xX ^ P{X) Reset map. 
In the above table, TX is the tangent bundle of X and P{X) is the power 
set of X. For more details on the definition of hybrid systems, see |14] . 
We can use the following representation for hybrid systems. 



where x G X is a vector of state variables and the form of f{x,X,q) is dic- 
tated by q, which represents the mode of operation of the hybrid dynamical 
system. The discrete state q is determined by the guard conditions that 
dictate transitions between modes (see Fig. [T]). The reset functions hi{x) 
are a part of the reset map R. In Eqn. [1] let A denote the vector of system 
parameters and x{0) the initial condition for the system. 

If the system parameters A in Eqn. [T] are uncertain (i.e., each Aj has 
an associated distribution) then one typically desires to quantify the time- 
varying moments (such as mean and variance) of x{t) (note that x(0) may 
also be uncertain). As mentioned earlier, although one can use Monte Carlo 
based sampling methods [l], they are plagued by slow convergence. In par- 
ticular, the mean is expected to converge as where N is the num- 
ber of samples. Quasi-Monte Carlo based sampling methods are expected 
to give a convergence rate of log'^(A^)/A^, where d is the dimensionality of 
the random space [24j . making these methods attractive for problems in 
low dimensions. Polynomial chaos based methods provide an alternative 
framework for uncertainty quantification with exponential convergence for 
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Figure 1: Schematic for hybrid (switching) systems. 

processes with finite variance, but they too suffer from the curse of dimen- 
sionahty [1]. In the next section we extend the polynomial chaos framework 
to hybrid systems. 

3 Polynomial chaos for hybrid dynamical systems 

Starting with a complete probability space T given by (0, J^, P), where is 
the sample space, J- is the cr-algebra on and P is a probability measure, 
let L2{T,X) denote the Hilbert space of square-integrable, F-measurable, 
X-valued random elements. Then one can, in general, define a polynomial 
chaos basis {Ha{X{tj))}, where A(a;) is a random vector and a = (ai, 02, • • • ) 
is a vector of non-negative indices. We denote the probability density func- 
tion of the random vector A by p(A). 

Generalized polynomial chaos (gPC) [25], provides a framework for rep- 
resenting second-order stochastic processes r S L2{T,X) for arbitrary dis- 
tributions of A by the following expansion: 

00 

r(A) = aaHaiX), (2) 

|a|=0 

where |a| = Oi is the sum of the indices of a and Ha{X) are orthogonal 
polynomials on T with respect to p(A), i.e. 

j^p{X)H^{\)Hp{\)dX = 6^p, (3) 

where 5ap is the Kronecker delta product. Depending on p{X) one can gen- 
erate an appropriate orthogonal basis for representing r(A). For example, 
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if p is Gaussian, then the appropriate polynomial chaos basis is the set of 
Hermite polynomials; if p is the uniform distribution, then the basis is the 
set of Legendre polynomials. For details on the correspondence between dis- 
tributions and polynomials see [11[26]. A framework to generate polynomials 
for arbitrary distributions has been developed in |25j . 

In practice, the expansion in Eqn. [2] is truncated at a particular order, 
say, p. One can then use Galerkin projections to obtain a set of differential 
equations for the coefficients aa in Eqn. [2] [3]. 

We now extend the standard polynomial chaos framework to hybrid dy- 
namical systems despite the presence of switching and state resets. To the 
best of our knowledge, it is the first attempt to develop tools for fast uncer- 
tainty quantification for this class of systems. 

3.1 Hybrid Polynomial Chaos 

Without loss of generality, consider the following two-mode hybrid dynam- 
ical system as representative of systems in which the different operating 
modes are associated with non-overlapping regions: 

\g{x^\) otherwise. 

Here one desires to quantify x{t] A), i.e., determine function of time t 

and parameters A. The system above has two modes of operation determined 
by its state. One can parameterize these modes in the following way: 

1 (x) - ifa;>0 

^ 1 otherwise 

1r,{x) = 1-1r,{x). (6) 

When (x) = 1 (corresponding to x > 0) the governing differential equa- 
tions are /(x,A), and when \ii^[x) = 1 (x < 0) the governing differential 
equations are g{x,\). Thus, one can rewrite Eqn. [4] as, 

X = lR^{x)f{x,\) + lii^{x)g{x,\). (7) 

This equation extends easily to k modes of operation by constructing indi- 
cator functions for each mode of operation , )•••)!/?* of the hybrid 
system. We now expand x{t] A) in the appropriate orthogonal polynomial 
chaos basis, 

p 

x{t-X)=Y,aa{t)H^{\). (8) 

|a|=0 
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Dropping the arguments of aa(t) and Ha{\) for simplicity and using the 
above relation with Eqn. [TJ one gets 



p 



p p 



a 



a\=0 



a\=0 \a\=0 



P P 



+ li?2( ^ aaHa)g{ ^ aaHa, A). 
la|=0 |o|=0 



By multiplying the above relation by /o(A)-fffc(A), integrating over T, and 
using orthogonality conditions, we get 



Note that Ri{t) = {A : Efa|=o'^"^" ^ ^2(i) = T - Thus 

by evaluating the two integrals one can evolve ak{t) for any index vector k. 
Note, however, that the regions of integration Ri and R2 are time-dependent 
quantities and must be evaluated at every instant in time. For a pictorial 
depiction of Ri and R2 , see Fig. [21 

3.2 Hybrid Polynomial Chaos and wavelet expansions 

Hybrid systems can display discontinuous behavior as a function of the un- 
certain parameters. In view of this, a smooth polynomial chaos expansion is 
expected to degrade as the discontinuities become more severe. In Ref. [22] 
the authors develop a wavelet-based Wiener-Haar expansion to treat bifur- 
cating (but smooth) dynamical systems with uncertain initial conditions 
that result in discontinuous behavior. In this section we adapt the Wiener- 
Haar expansion to hybrid dynamical systems. 

In \22\ [23], output variables are expanded in terms of Wiener-Haar 
wavelets expressed as functions of the Cumulative Distribution Function 
(CDF) of the uncertain parameters. For simplicity, consider the univariate 
case. Here we denote the CDF of the uncertain parameter A as u{X) and 
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Figure 2: Pictorial representation of the regions of integration for the two 
integrals in Eqn. [9l In this example, the regions are not simply connected. 



expand the state vector X a£0 



P 2^-1 

x{t; A) = xo{t) + ^ ^ Xjk{t)ipjk{u{X)), 

j=0 k=0 



(10) 



where, 

ipjk{u) = 2^/^(2^'" - k), (with J = 0, 1, . . . and /e = 0, . . . , 2^' - 1) 
is a family of Haar wavelets [21], defined in terms of the mother wavelet: 



1 0<n<l/2 
-1 l/2<u<l 
otherwise. 



The index j determines the scale of the wavelet and k its displacement. Note 
that {tpjk} is a family of orthonormal functions on the interval [0, 1] with 
respect to the uniform density. This makes the family {ipjk°u} automatically 
orthonormal with respect to the probability density of A: 

SjiSkm = j i^jk{u)i)im{u)du = j {ipjk o u){X) o u){\) p{X)dX. 



^For the multivariate case, see Ref. 



7 



Additionally, all ipjkS are orthogonal to the constant function on [0, 1], which 
implies that the mean of x is given by the first term in the expansion 



Xo{t) = J x{t;X)p{X)d\ 

and that the variance is 

We now use this expansion on a switching oscillator example: 
x + cx + x + X = ifx>0 

x + cx + x-\ = ifx<0, (11) 
which can be rewritten as, 
X = y 

y = -cy -X - A1r,(x) + \1r^{x). 
The expansion in this case is 

P 2J-1 

x{t]X) = Xo(t) H ^3k{t)'^jk{u{\)) 

j=0 fc=0 
P 2^-1 

y{i] A) = 2/o(t) + X] X] y3k{t)i^jk{u{X)). 

j=0 k=0 

Projecting these equations onto the basis functions yields, 

xq = yo 

m = -cyo -xq- X{u)lR^{x{u))du + / \{u)lR2{x{u))du 
Jo Jo 

•^jk Vjk 

iljk = -cyjk - Xjk - / X{u)lR^{x{u))'il)jk{u)du + / \{u)lR^{x{u))tl)jk{u)du 
Jo Jo 

Note that to compute \{u) we invert the CDF u{\). 

To compute the integrals needed to evolve these equations numerically, 
we take advantage of the fact that Haar wavelets are piecewise constant. 
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Namely, for a given truncation order P, if we divide [0, 1] into 2'^'^^ equal 
subintervals, both ■0^^ (with j < P) and the truncated expansion for x (and 
therefore the indicator functions) are constant in each subinterval. This 
implies that in each subinterval we only need to calculate the integral of 
A(ti), which is known a priori. For the case of a Gaussian A ~ iV(/x,(T^), we 
have 

X{u) = fi + aV2err^{2u - 1), 

which has a primitive, 

J X{u)du = fiu — o—-j= exp | — [erf~"'^(2n — 1)]^} . 

Therefore the contribution to y^fe of the integrals in each subinterval / = 
0, . . . , 2^+^ - 1 is either zero or ±2-'/^ times the precomputed value 

(/+l)/2^'+l 

\(u)du. 

Section 13.41 presents the results obtained with the Wiener-Haar wavelet ex- 
pansion in Eq. [10] for hybrid dynamical systems. 



3.3 Modeling state resets 

A significant challenge that hybrid dynamical systems present is the possi- 
bility of state resets jl4] . When a hybrid system switches from one mode 
to another, the state of the system can, in general, be reset discontinuously. 
For example, in the case of the bouncing ball [13], the velocity of the ball 
changes discontinuously after every impact. When the hybrid system tran- 
sitions from mode q to q' the state resets are typically represented as, 

x+ = h{x-), (12) 

where x~ and x"*" are the states of the system before and after the reset. Such 
discontinuities cannot be easily accommodated within the hybrid polynomial 
chaos framework as described in the previous sections. 

To circumvent this problem, one can construct a boundary layer in the 
vicinity of the guard condition. We also introduce a dummy vector (z) that 
tracks the state x outside the boundary layer and is set to x~ within the 
boundary layer. Note that we assume separability of the states, i.e. the 
guard conditions (which determine the switching between modes of opera- 
tion) can be written independently of the state reset conditions in Eqn. [T2l 
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In other words, the states that determine the guard conditions do not par- 
ticipate in the state reset. Let the reset condition be in terms of vector 
X in Eqn. [T2l and the guard conditions be in terms of vector y (given by 
{y '■ diy) = 0}). Note that, [x y]^ represents the entire state vector with 
the following governing equation, 

X = fi{x,y) 

y = f2{x,y). (13) 

We now construct a boundary layer around the guard condition for vector 
y as follows: 

if:\9{y)\>e 

(14) 

otherwise. 

The above dynamical system is constructed such that x~ evolves to 
h{x~) in At ~ e, where e is a small parameter. 

By replacing each reset condition with an equation of the form given by 
Eqn. [ni one obtains a new dynamical system without resets that approxi- 
mates the original. On this new dynamical system one can use the expan- 
sion from Sec. 13. II and evolve it using Eqn. [9l In other words, the framework 
generalizing polynomial chaos to hybrid systems can be augmented using 
Eqn. [T3] to include state resets. 

To illustrate the procedure presented above we turn to the classic bounc- 
ing ball example [15] : we consider the dynamics of a ball bouncing on a floor 
with coefficient of restitution 7 < 1 under the action of gravity of uncertain 
magnitude {^jL{g) = 9.8 m/s^ and a{g) = 0.2 m/s^). Thus, every time the 
ball makes contact with the floor the velocity v~ is reset to a new value 
given by f"*" = —^v~. The guard condition for resetting the velocity is given 
by y{t) = (where y{t) is the height of the ball above the floor at time t). 
The equations for the bouncing ball are given by, 

y = V, 

V = -g, (15) 



f h{x.y)\ 
f2{x,y) 
\{x - z)/ej 



V 



f[h{z)-x]/e 
^f2{x,y) 
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Figure 3: a) Monte Carlo simulation of the bouncing ball system with uncer- 
tain gravitational acceleration, b) Nominal bouncing ball trajectory com- 
pared with the mean trajectory obtained through the wavelet-based hybrid 
UQ approach with the boundary layer approximation (e = 0.01). 
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with the reset condition at y = 0: v'^ = —jv 
boundary layer approximation in Eqn. [T3]we get, 



Thus, if one uses the 



r / 



if : |y| > e 



\{v-z)/e, 



(16) 




otherwise. 



We now use the hybrid polynomial chaos expansion in Eqn. [9] along with 
the Wiener-Haar wavelet basis functions. The Monte Carlo simulations on 



the bouncing ball are shown in Fig. 3(a) , The average or nominal trajec- 



tory is also shown. In Fig. 3(b) we compare the nominal trajectory (mean 



trajectory from Monte Carlo) with the mean predicted using the boundary 
layer expansion with e = 0.01. As shown, the boundary layer accurately 
approximates the mean over multiple state resets events (in this case, each 
impact with the floor). 



3.4 Results 

To demonstrate the hybrid polynomial chaos approach on hybrid dynami- 
cal systems we consider the simple yet challenging example of a switching 
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oscillator given by Eqn. [TTJ 

x + cx + x + X = ifx>0 
x + cx + x — X = ifx<0. 

The value of c is deterministic and equal to 0.5. Here we consider three 
cases with A normally distributed with: /x(A) = — 10 and cr{X) = 2 (case 1), 
/x(A) = 10 and a{X) = 2 (case 2), and //(A) = and a{X) = 1 (case 3). In 
all cases we assume that the initial conditions are deterministic and given 
by [x{0),x{0)] = [10-2,1.0]. 

3.4.1 Case 1: ^(A) = -10, a{X) = 2 

Let us start with the case when /-f(A) = — 10 and o"(A) = 2 in Eqn. [TTJ A rep- 



resentative trajectory for the dynamics of the system is shown in Fig. 4(a 



The corresponding histogram for x(20.0; A) is shown in Fig. 4(b) Most 
importantly, one desires to compute the mean and variance of x{t; A) as a 
function of time. In the system given by Eqn. \TT\ we expand x{t; A) using 
Eqn. [8] and perform a Galerkin projection as shown in Eqn. [9l One then 
gets a system of equations for the coefficients of expansion in Eqn. [8l These 
coefficients, once computed, can be used to calculate the moments of the 
distribution of x{t;X). 

We compare the results obtained from hybrid polynomial chaos with 
those obtained using Monte Carlo and Quasi-Monte Carlo based methods. 
In particular, we use a Weyl sequence [27] along with inverse transform 
sampling [28j to generate the Quasi-Monte Carlo samples. We find that the 
results (in the first two moments) from 5000 samples of Monte Carlo, 3000 
samples of Quasi-Monte Carlo and the Wiener-Haar hybrid PC expansion 



with P = 3 are visually indistinguishable (see Figs. 5(a) and 5(b) ). Treating 
5000 Monte Carlo samples as baseline, we find that the hybrid PC expansion 
has a maximum error of 5 x 10"^ in the prediction of /i(x(t; A)). 

3.4.2 Case 2: ^(A) = 10, cj(A) = 2 

The case of /i(A) = 10 and cr(A) = 2 is significantly more challenging. A 



representative trajectory of the system (for A = 10) is shown in Fig. 6(a 
When A > the system switches back and forth between modes. The reason 
for this is that when the system is in the right half-plane, the equilibrium of 
the system is in the left half-plane and vice versa. A histogram for x(3.0; A) 
is depicted in Fig. |6(b)[ 
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Figure 5: Comparison of a) predicted mean and b) predicted variance of 
x{t; A) for various UQ methods for case 1. 
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We again compare hybrid Wiener-Haar polynomial chaos to Monte Carlo 
sampling in Fig. [71 We find that hybrid Wiener-Haar polynomial chaos 
{p = 5) accurately computes the mean and the variance of the distribution 
of X. The maximum absolute error of hybrid polynomial chaos in mean is 
fj,{x{t; X) = 1.8 X 10~^ and variance is 7 x 10~^. Note that expansions in 
terms of standard basis functions such as Hermite and Legendre polynomi- 
als are unable to compute the moments of x{t; A) beyond a threshold time 
that depends weakly on the order of expansion p (see Fig. [8]) . The solution 
in this case is particularly challenging because it becomes more oscillatory 
in terms of A at t increases (Fig. [9]). The Wiener-Haar basis functions are 
naturally oscillatory and hence more accurate than Hermite polynomials 
in capturing the solution x{t;X). Note that, for large time simulations the 
Wiener-Haar expansions will also fail since the solution will eventually be- 
come too oscillatory for the order of expansion. This problem is well known 
in the polynomial chaos literature [29]. 



3.4.3 Case 3: ^(A) = 0, cj(A) = 1.0 



We also consider the case of /x(A) = with cr(A) = 1.0. This case is partic- 
ularly challenging because there is a concentration of probability of x{t; A) 
as shown in Fig. 10(b) , The reason for this is as follows: when /x(A) = 0, 
the nominal trajectory converges to and so do all trajectories with A > 0. 
Indeed, for trajectories with A > the equilibrium lies in the opposite half- 
plane with respect to the current state. This gives rise to decaying switching 
trajectories, as case 2 in Fig. 6(a) Note that for A < 0, the trajectories are 
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Figure 7: Comparison of a) predicted mean and b) predicted variance of 
x{t; A) by various UQ methods for case 2. 
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Figure 8: a) Mean and b) Variance predicted by standard (Hermite) poly- 
nomial chaos basis functions for case 2. 
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Figure 9: Case 2: x{t; A) becomes more oscillatory in A as t increases. The 
bottom plot shows the distribution for A. 
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Figure 10: a) Nominal trajectory for case 3. b) Histogram of x(20;A) for 
case 3. 



similar to the ones in case 1 (Fig. 4(a) ). 

The Wiener-Haar basis functions along with the hybrid polynomial chaos 
approach accurately capture the moments of the distribution for x{t; A) (see 
Fig. Ilip. In fact, an expansion to just P = 3 captures the first two moments. 
The step-function nature of the Wiener-Haar basis allows it to perform well 
in this scenario. Standard basis functions like Hermite polynomials are com- 
pletely incapable of accurately capturing the moments of the distribution for 
x(t; A) shown in Fig. |6(b)[ 



4 Transport theory approach for uncertainty quan- 
tification in hybrid systems 

In this section we present a qualitatively different approach to UQ in hy- 
brid systems based on transport equations. We write an advection equation 
for the probability density of the state and expand this equation in an ap- 
propriate basis, as is done in polynomial chaos. The resulting equation is 
equivalent to the Fokker-Planck equation |30] in the absence of a diffusion 
term. Though significant effort has been put into computing solutions for 
the Fokker-Planck equation in various applications [30\ [31] . our setting is 
particularly challenging due to the switching dynamics of hybrid systems. 
We note that advection equations for probability distribution functions have 
been used to propagate uncertainty through heterogeneous porous media 
with uncertain properties |32j and for hyperbolic conservation laws with 
noise |33] . Recently, similar methods have been extended to cumulative 
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Figure 11: Comparison of a) predicted mean and b) predicted variance of 
x{t; X) for various UQ methods for case 3. 



distribution functions in hyperbolic conservation laws |34j . 

The polynomial chaos expansion in this setting yields a system of hyper- 
bolic partial differential equations for the coefficients of the expansion, which 
are then solved by integrating along characteristics. The hybrid nature of 
the original system is reflected in that the characteristics exhibit switching. 
Even though we only consider systems without resets, we can use the results 
of Sec. IS.SI to treat systems with resets. 

As in SecO let us consider a hybrid system without resets and uncertain 
parameters A with guard conditions independent of aH 

X = fi{x, X) when Gj(x) is true. 

The system has uncertain initial conditions described by the probability 
density pxo{x) and the uncertain parameters A follow px{X). 

We describe the system by the time evolution of the distribution function 
p{x,X;t), which has initial condition 

p{x,X;0) = Pxo{x)px{X) (initial uncertainties are independent) 

and normalization 

J p{x, X;t)dxdX = 1. 

Note that, for all time, we have 

px{X) = j pix,X;t)dx. (17) 

^Note that this embodies the constraint of having no overlap in the domains for different 
modes. 
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Our goal is to compute the evolution of the density in x (the marginal 
distribution) : 

Px{x,t) = j p{x,X;t)dX. 

However, without introducing assumptions on px, the equation for is 
not closed. We therefore focus on computing the evolution of p directly 
through an expansion. From this evolution, p^ can then be calculated at 
every instant. 

4.1 Equation for p 

Let us define the sets Si = {x : Gi{x) is true} and the indicator functions 



1 \i X ^ Si 
if not. 



With this notation, and because A is constant along a trajectory, we have 

|^ + V-(p/) = VA (18) 

where f{x,X) = lj(x)/i(x, A) and the gradient operator acts only on x 
and not on A. 



4.2 Boundary conditions at the interfaces 

The discontinuity in the equation implies that mass may accumulate at 
the boundaries between zones where different guard conditions are valid. 
Integrating on a cylinder that crosses one such boundary we obtain the 
matching condition 

0(7 

— + Vs • (fj/) = Pifi • Uik - Pkfk • nik, (19) 

where o" is a surface probability density between regions Si and Sk, Uik is 
the surface normal from Si to Sk, Vs is the divergence in the space tangent 
to the surface, and / is the flow at the surface^ This may lead to a cascade, 
with probability condensing into progressively lower dimensional structures: 
where two hypersurfaces meet (the boundary between three guard condi- 
tions) the same scenario repeats, until we have mass accumulating at points. 

^Whether / is fi or on the surface will depend on how the guard conditions are 
expressed. 



19 



To solve Eqn. [18] the initial condition must include initial values for a and 
the probability density on any lower dimensional structure where mass may 
accumulate. 

The discontinuity of / does not necessarily imply accumulation. In fact, 
at an interface we have several options: 

1. Both fi ■ hik and fk ■ hik are nonzero and have the same sign. In this 
case, there is no accumulation. If we assume that p has a singularity 
at the interface, the flow will move the singularity away from it. 

2. fi ■ hik > and fk ■ hik ^ 0: accumulation occurs. 

3. fi • hik > and fk • hik < 0: accumulation occurs. 

4. fi ■ hik = fk ■ nik = 0: no accumulation. 

5. fi ■ hik ^ and fk ■ hik ^ 0: no accumulation. 

If we are in the case without accumulation and without initial concentration 
of density in lower dimensional structures, then a = and Eqn. 1191 becomes 



satisfies the conditions for no accumulation at the interface where the ODE 
is discontinuous. 

Proof. Without loss of generality we can focus on just two regions Si and 5^ 
and rewrite the problem as the first-order ODE 



To find the normal h^k we consider a function 6(x) = b{X,Y) = B{X) 
that is positive in Sk and negative Si so that the interface is given by the 
locus of b[x) = 0. The gradient of this function is proportional to the normal: 



Pifi ■ flik = Pkfk ■ flik- 

Theorem 1. Any second order ODE of the form 

X = Fi{X,X, A) when Gi{X) is true 



(20) 




when Gi{X) is true. 




and therefore 



flik ■ f 



l|Vx-B|| 



which is continuous at the interface. 



□ 
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4.3 Expansion of the equation for p 

At every point x we expand the distribution in A: 

p{x, X;t) = Y^ ak{x, t)w{X)iPkW, (21) 

k 

where {V'fc} forms an orthogonal basis with respect to w: 
j '4)i{\)ilJk{X)w{\)d\ = Wk5ik- 

We keep Wk to allow for a non-normalized weight function w{\). We replace 
the expansion in Eqn. [21] into Eqn. [18] and project onto tpi to obtain a set 
of partial differential equations for the coefficients ai{x,t): 

+ ■ ^ / ^(^)V'.(A)V'fc(A)/(x, X)d\ = 0. (22) 

k 

Since this equation is local in x there is no question as to which f{x, A) must 
be used at any given point. 

4.4 Example: switching oscillator 

Here we revisit the switching oscillator system 

{—X — jx — X if X > 
—X — 7x + A otherwise, 

which can be expressed as the 2-D system 

© ^ (t) 

where fx = y and 

I —X — 7y — A if a; > 



"^^ 1 — X — 7y + A otherwise. 
The transition points are located at x = and therefore 

f ■ n = f ■ = fx = y 

which is continuous and therefore has the same sign on both sides. Therefore 
there is no mass accumulation at the interface for this system. This is a 
special case of theorem [T] 
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4.4.1 Case 3 revisited 

To connect with the example presented in case 3 (/x = 0, cr = 1) we choose 
w{X) = e~^^ 1"^ and i/'fc's as the probabilist's Hermite polynomials i/fc, with 
the following properties 



Hi{X)Hk{\)w{X)dX = kW27r6ik 

Hk+i{X) = XHkiX) - H'j,{X) 

Calculating the terms in Eqn. [22l 
wHiHkfxdX = yi\V2TT6ik 

wHiHkfydX = -{x + jy) ilV27r6ik =F i!\/27r(5i,fc+i + k6i^k-i) 



I 



where the upper sign (— ) is for x > and the lower sign (+) is for x < 0. 
Substituting into Eqn. [22] we obtain the equations 

= dtao + dxiyao) + dy [-{x + 7y)ao T ai] ^^3) 
= dtai + dxiyai) + dy [-{x + ^y)ai =F (i + =F Oj-i] («>!)• 

Note that, instead of using the Hermite polynomials, one can use the Haar 
wavelet expansion [22] to represent the solution as was done in previous 
sections. We plan to present this calculation in future work. 

Theorem 2. The system of PDEs for the switching oscillator is hyperbolic. 

Proof. Consider a system of PDEs of the form 

dtu + ^ AuduU = B, 

where u(xi, . . . , x„,, t) G M™ and the A^, are m x m matrices. The system is 
hyperbolic if for any £M the linear combination A = Oi^Ay has real 
eigenvalues. 
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For the switching oscillator the system of PDEs can be written as 



dt 




(y .... 




/ao\ 


y 




ai 


V 




vl 










/3 



+ 



Tl /? =F3 



v ■■ 









/ao\ 




ai 




ai 











where /3 = — (x + 7y). Thus, any combination of the matrices Ay is going 
to be of the tridiagonal form 

0\ 



A 



a 


b 







b 


a 


26 








b 


a 


36 







6 


a 



■I 



This tridiagonal non-symmetric matrix is similar to a tridiagonal symmetric 
matrix with a diagonal similarity matrix: S = DAD~^, where 



D 



\ 







V 







•••/ 



A is similar to S, a symmetric and real matrix, and therefore A has real 
eigenvalues. 

□ 

The issue of hyperbolicity is discussed in more depth in [35]. To solve 
the hyperbolic system from Eqn. [23] we write it in the form 



dta + ydxtt — (x + ^y)dya =F AdyU = ja, 
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where A is the tridiagonal matrix 

/O 1 



A 



10 2 

1 3 





1 



V 



(24) 



•7 



We diagonahze A = PAP ^ and define b = P to obtain the set of 
uncoupled hyperbohc PDEs 



dfbi + ydxbi + [-(x + 7y) =F Aj] dybi = 'jbi 



(25) 



where Aj is the i-th eigenvalue. We will now prove that when the expansion 
is truncated up to a„_i, i.e., A is truncated to a nxn matrix, the eigenvalues 
A are the zeros of H^- 

Theorem 3. The eigenvalues of An, the nxn truncated version of the 
matrix in Eqn. \24[ are the zeros of the n-th order probabilist's Hermite 
polynomial Hn- 

Proof. We proceed by induction to prove that det{An — XI) = i?„(— A), 
which will then, by the symmetry of Hn, prove our result. 

Let Bn = An — XI. Indeed, det Bi = —X and det B2 = - 1. For the 
general case, 

/ o\ 



B^ 



n+l 



Br, 





n 



1 



-A/ 



Therefore, 

detS„+i = -AdetS„ + (-l)"n(-l)"-MetB„_i = -A det B„ - n det 
which is the recurrence relation satisfied by Hn{—X). 



□ 



The characteristic curves of Eqn. [25] are given by 

x = y 

y = -X - 71/ =F Ai 
hi = -fbi. 
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In other words, the characteristics are damped oscillators where the equilib- 
rium position is given by the eigenvalues of ^A. The exponential growth of 
bi along a trajectory is due to the contraction in phase space produced by 
the dissipation 7. 

4.4.2 Results 

We now show results obtained by using transport theory approach on case 3 
(/U = 0, fj = 1) for the switching oscillator (Eqn. [TT|) . 

In Fig. [T2I we show a series of probability distribution snapshots for 
Monte Carlo (5000 samples) and the transport operator method (for case 3), 
gridded in the (x, y) plane. Note that we set y = x, as defined in the first part 
of the paper. The Monte Carlo color map snapshots show that the distri- 
bution lies on a one-dimensional manifold (in two dimensional space). The 
one-dimensional nature of the distribution arises because we chose deter- 
ministic initial conditions. As discussed previously, all trajectories in case 3 
with A > converge to the origin, resulting in a jump in the cumulative 
distribution function (CDF). 

For A < 0, however, the trajectories converge asymptotically to either +A 
or —A. This convergence to +A or —A is highly dependent on the individual 
trajectory and results in a fragmentation of the output distribution. Close 
to convergence {t = 18), very few trajectories converge to a point in the 
range 0.2 < x < 0.4, as seen in the fiat region of the CDF in Fig. [T3j The 
transport-based method captures the singularity at the origin accurately, 
but is unable to accurately capture the fragmentation. This is because 
the method samples the distribution sparsely (determined by the order of 
expansion), resulting in UQ acceleration. However, this sparsity makes the 
method miss such fine details. On using a high order of expansion (n = 70), 
some samples partially capture the structure around x = 0. Note that a 
much lower order expansion accurately captures the jump at the origin and 
the asymptotic (x > 1) shape of the CDF. 

As in Fig. [TT| we compare Monte Carlo (5000 samples) with the trans- 
port theory approach (orders 15 and 70 expansion) in Fig. [T4l The h^o 
(maximum) errors for n = 15 are 9.56 x 10~^ (mean) and 7.95 x 10~^ (vari- 
ance). For n = 70 we get Lqo errors of 5.28 x 10"^ and 4.92 x 10"'^ for mean 
and variance respectively. As shown in Fig. [TTl the method performs rea- 
sonably well, however, the results are not nearly as good as those obtained 
using hybrid polynomial chaos with the Wiener-Haar wavelet expansion in 
section [Xn However, with a better choice of basis functions, one does expect 
better results. The transport operator theory is attractive as it appears to 
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Monte Carlo 5,000 samples Transport theory approach n = 70 CDF 




XXX 



Figure 12: Snapshots at t = 2, 4, 6, 8, 10 of the gridding from 5000 sam- 
ples of Monte Carlo and an order 70 expansion using the transport theory 
based method. The right column compares the one-dimensional cumulative 
distribution function for the corresponding times. 
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Figure 13: Comparison at t = 18 of the cumulative distribution function 
from 5000 samples of Monte Carlo (solid) and an order 70 expansion using 
the transport theory based method (dashed). 



be more versatile. In general, the transport operator approach is applicable 
to hybrid dynamical systems with overlapping modes of operation (by con- 
structing multiple PDEs for the overlapping mode). In contrast, the hybrid 
polynomial chaos method suffers from the disadvantage of being inapplicable 
to such systems. 



5 Conclusions 



As the modeling of hybrid dynamical systems becomes increasingly impor- 
tant for modern day engineering applications such as electrical and biologi- 
cal networks, air traffic systems, communication networks, etc., quantifying 
uncertainty in these systems is going to become a central concern. Since 
uncertainty quantification allows one to compute moments of output distri- 
butions in the presence of parametric uncertainty, these techniques will be 
used to aid decisions related to robust system design and performance. 

In this work, we have made the first attempts to develop fast uncertainty 
quantification methods targeted for hybrid dynamical systems. In particu- 
lar, we extended polynomial chaos methods, a popular technique for prop- 
agating uncertainty through smooth systems, to hybrid dynamical systems. 
We also developed methods to handle state resets within the polynomial 
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Figure 14: Comparison of a) predicted mean and b) predicted variance of 
x{t; A) for case 3. The curves show a 5000-sample Monte Carlo run, and 15 
&: 70 term expansions using the transport theory approach. 

chaos framework by using boundary layer approximations. We then apphed 
this new approach to perform uncertainty quantification on switching har- 
monic oscillators and the bouncing ball examples. We also demonstrated 
the efficacy of using Wiener-Haar expansions [9l [22| [23] with our hybrid 
polynomial chaos approach for quantifying uncertainty in hybrid systems 
that give rise to multi-modal distributions or become increasingly oscilla- 
tory in time. Finally, we showed how a transport theory based approach 
can capture naturally-emerging discontinuities in the distribution. Future 
efforts involve providing rigorous error bounds for Wiener-Haar expansions 
in the hybrid polynomial chaos setting with boundary layer expansions. We 
are also extending our hybrid polynomial chaos approach to networks of hy- 
brid dynamical systems using our recent work on propagating uncertainty 
through complex networks [36]. We also intend to extend the transport 
operator based UQ method to hybrid systems with overlapping modes of 
operation. 
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